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ALTARELLI-PARISI EVOLUTION EQUATIONS* 
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Abstract: We discuss a new method to solve in a semianalytical way the Dok- 
shitzcr-Gribov-Lipatov-Altarelli-Parisi evolution equations at NLO order in the ir-space. The 
method allows to construct an evolution operator expressed in form of a rapidly convergent series of 
matrices, depending only on the splitting functions. This operator, acting on a generic initial distri- 
bution, provides a very accurate solution in a short computer time (only a few hundredth of second). 
As an example, we apply the method, useful to solve a wide class of systems of integrodifferential 
equations, to the polarized parton distributions. 
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1. Introduction 

The scaling violation of nucleon structure func- 
tions is described in terms of Dokshitzer-Gribov- 
Lipatov-Altarelli-Parisi (DGLAP) evolution equa- 
tions [1]. The DGLAP integrodifferential equa- 
tions describe the Q 2 dependence of the struc- 
ture functions, which are related, via the oper- 
ator product expansion, to the parton distribu- 
tions for which the DGLAP equations are usu- 
ally written down. In this framework, the analy- 
sis of the experimental data, is performed fixing 
at some <5o the structure functions by assuming 
the parton distributions and computing the con- 
volution with the coefficient functions, which can 
be evaluated in perturbation theory. The com- 
parison with experimental data, which are dis- 
tributed at different values of Q 2 , goes through 
the solution of DGLAP equations for the parton 
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distributions; thus a reliable and fast algorithm 
to solve these equations is welcome. 

In literature there are essentially three differ- 
ent approaches to solve the DGLAP equations. 
The first one is based on the Laguerre polynomi- 
als expansion [2]. This technique is quite accu- 
rate up to x- values not smaller than x « 10 -3 ; on 
the contrary, below x the convergence of the ex- 
pansion slows down [3, 4]. Given that experimen- 
tal data are already available down to about x, 
for the polarized case, and down to 1CP 5 , for the 
unpolarized case, this method results no longer 
practical. 

An alternative approach takes advantage of the 
fact that the moments of the convolutions ap- 
pearing in the equations factorize in such a way 
that the analytical solution, in the momentum 
space, can be written down [5]. However, also in 
the most favorable case in which the analytical 
expressions of the moments of the initial condi- 
tions arc known, the numerical Mellin inversion is 
relatively CPU time consuming (see [6]). More- 
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over, as discussed in [7], since x variable is related 
to the invariant energy W 2 of the virtual photon- 
hadron scattering process by W 2 = (1 — x)/x, 
x — > is the infinite energy limit and thus can 
never experimentally be reached. As a conse- 
quence of this all moments are plagued by an 
a priori infinite uncertainty, which can be re- 
duced by means of assumptions implying that 
any use of the evolution equations for moments 
is model dependent. The more simple solution to 
this problem is to solve DGLAP equations in the 
x-space. In this framework, besides the Laguerre 
method, another strategy, the so called "brute 
force" method [8], represents a good candidate. 
It is fundamentally a finite-differences method of 
solution which reaches a good precision in the 
small x-region [9, 10] but requires a rather large 
amount of computer running time. 

Here we discuss a semianalytical method, in 
the x— space, to solve DGLAP equations [11]. 
It consists in constructing an evolution opera- 
tor which, depending only on the splitting func- 
tions, can be worked out once for all. In this re- 
spect our strategy is similar to the one in [2, 4]. 
Our method to perform the convolutions instead, 
takes advantage of an x discretization (compara- 
ble to the one in [9, 10]) which allows us to rep- 
resent the evolution operator as a matrix. Thus 
the procedure to construct the solution reduces 
merely to a multiplication between the evolution 
matrix and an initial vector, and can be done in 
an extremely short computer time with the re- 
quired accuracy. This is particularly appealing 
in the analysis of the experimental data on nu- 
cleoli structure functions which requires a large 
number of parton evolutions. 
In the next section we discuss the (formal) ana- 
lytical solution of the DGLAP equations; in the 
third one the algorithm to perform the x-integration 
is presented. The last two sections are devoted to 
analyze the numerical results relative to the evo- 
lution of polarized parton distributions, to study 
the yield of our method in comparison with oth- 
ers and to conclude. 

2. The Evolution Operator 

The DGLAP equation, up to Ncxt-to-Lcading- 
Order (NLO) corrections, for the Non-Singlet dis- 



tribution is* 

— Aq NS (x,t) = 



AP { °l(x) + a(t)AR NS {xf) ® Aq NS (x, tj)2.1) 

while for the Singlet and Gluon distributions we 
have: 
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where 



AR ij {x)^AP^\x)-^-AP§\x). (2.3) 
In these equations the symbol <S> stands for 

f{x)®g{x)= f ' *Lf (*) g(y) , (2.4) 

Jx y \yJ 

and 

f(x) = xf(x). (2.5) 

Instead of Q 2 , we have used the variable t defined 
by 

■a s (Q 2 ) 



2 , 
In 



(2.6) 



where a s is strong running coupling constant cor- 
rected at NLO: 



ocs{Q 2 ) 



4ir 



Po ln(Q 2 /A 2 ) 
and so in Eqs. (2.1)-(2.2) 



A ln(ln{Q 2 /A 2 )) 
(3 2 ln{Q 2 /K 2 ) 



a(t)^^Ex P {-^t 



(2.7) 



(2.8) 



The explicit expressions for /3 , /3i as well as for 
the Splitting Functions APij(x) can be found in 
[12]. 

The equations (2.1) and (2.2) can be written 
in the following general form: 



|f(i) = n(t)sf(t) 



(2.9) 



*In the following we limit ourselves to discuss the 
polarized parton distributions. The application of our 
method to the unpolarized case is straightforward. 
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where f(t) indicates the "vector of components 
f(x,t)" and fi(t) a linear operator acting as: 

[«(t) Qf{t)] x = [ dy u(x,y,t) f(y,t). 

(2.10) 

Note that, in the Singlct-Gluon case, f(t) be- 
comes a doublet of vectors and fi(t) a 2 x 2 ma- 
trix of operators. 

Due to the logarithmic dependence of t on 
Q 2 , the range of values of physical interest for 
t — to (t is the starting values of t, where the 
parton distributions are assumed known) is small 
enough to expect that the Taylor's series of the 
solution f(f) converges rapidly. On the other 
hand, by deriving repeatedly the Eq. (2.9) we 
can write: 



gk 



M«Of(t ), (2.11) 



t=t 



where the operators M( fc ) can be obtained recur- 
sively: 

M<°> = I 

m (1) = n Q 

M (2) = JI^ + OoOMW 

m (3) = n 2) + 2 n 1} m (1) + n © m (2) 



fe-1 



M (fc) = 



(2.12) 



i=0 



The indicates the i— th term of the fc— th row 
of Tartaglia triangle and 



ft = n(to) 



5 fe 



n(t) 



*=*0 



(2.13) 



Then the solution can be written as: 

f(t) : 



.f(/„) 



fc! 

\fe=0 

T(t - t ) f (to) 



(2.14) 



with T(t — t ) the Evolution Operator. As we 
will point out in section 3 the series in Eq. (2.14) 
converge quickly enough to obtain a very good 
approximation retaining only a first few terms. 
It is worth to note that if the operator O(t) can 
be written as h(t) f2' (with h(t) a numerical 



function) it is easy to show that the series in 
Eq. (2.14) reduces to: 



f (t) = Exp ■ 



{[/: 



h(r)dT 



O'lof(to). (2.15) 



This is the case of DGLAP equation at Lead- 
ing Order (LO) approximation. Nevertheless, in 
Eqs. (2.1)-(2.2), where NLO corrections arc in- 
cluded, we have fi(t) = fli + a(t)fi 2 , with fii 
and $^2 non-commuting operators. As a conse- 
quence the series in Eq. (2.14) cannot be summed 
and it is not possible write the solution in a closed 
form. 



3. The ^-Integration 

The integrals in Eq. (2.4) are evaluated with a 
method that generalizes the one proposed in Rcf. 
[10]. The method consists to treat exactly the 
"bad" behaviour of the kernel w(x, y, t) in Eq. (2.10) 
and approximate the "smooth" function f(y,t). 
In particular, we construct a M + 1 points grid 
(xo > 0, xi, Xm-i, xm = 1) in the interval 
]0, 1] and approximate f(x) in each interval [xk, Xk+i] 
by the cubic which fits the four point f(x{), with 
i = k-l,k,k+\,k + 2: 

4 

f(x) w ^2 a[ k) x l ^ 1 (x) Vx G [x k , x k +i] ■ 
i=i 

(3.1) 

The general structure of the Polarized Splitting 
Functions which appear in Eqs. (2.1)-(2.2) is t : 



A{x) 



2) 



AP(x) = — M- + B(x) + 6(1 - x)C , (3. 

(1 - x) + 



and therefore the "i component" of the convolu- 
tion is: 



fl dyA(x i /y)f(y)-A(l)f(x i ) 



y 



y- %i 



+ 



/>(?M + 

(C + A(l)ln(l-x i ))f(x i ) . 



(3.3) 



tThe same structure, however, holds for unpolarizcd 
and transversely polarized splitting functions. 
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Substituting Eq. (3.1) in Eq. (3.3) we obtain Vi £ 



{0, M - 1} ( Y!k=M = is understood) 

m 

&p(xi) ® f{xi) = y,4 ] (Pi+pii) + 
i=i 

M-l m 

E E a ; (fe) (^+^) + 

fe=i+i ;=i 

(C + ^(l)/n(l - x,) - ,4(1)*') /(a*) ; (3.4) 
then we have: 

M 

AP(x i )®f(x i ) = J2"ikf(xk), (3.5) 



fe=0 



where u is the matrix of the coefficients of f(xk)- 
The analytical expressions for the matrices f3, 7, 
p and (T can be found in [11]. 

Therefore the Eqs. (2.1)-(2.2) became *: 

3 A - / \ 

—Aq NS {x ll t) = 

M 

E(4 0) AfS + ^K ( fc 1) WS ) ^(^,0(3.6) 

fe=0 



^ ^ A3 (a*,*) 



M 



= E 



, ,(0) 9 g 95 
, ,(0) 39 (0) g S 



( , ,(1) 11 , ,(1) 9ff ' 
+«(*) ^) ,9 > S9 



(.3.7) 



We solve these equations by means of the method 
shown in section 2: the operator Sl(t) and then 
the M( fc ) became now numerical matrices, and 
the symbol stands for the usual rows by columns 
product. We would stress the fact that the ma- 
trices M( fe ) depend only on the points Xi and so 
they can be numerically evaluated once for all. 

4. Numerical Analysis 

The convergence of our algorithm is controlled 
by two parameters: the order n of the truncated 
series 

T (")(*-to) = E^^MW , (4.1) 



tNote that a/ 1 ) matrices correspond to the convolu- 
tions of the parton distributions with AR (cf Eq. (2.3)). 



which define the evolution operator, and the num- 
ber M of the points of x— integration. 

To test the accuracy of our method we evolve 
the Gchrmann and Stirling polarized singlet-gluon 
initial distributions (cf [13]) from Ql = 4 GeV 2 
(to = 0) to Q 2 = 200 GeV 2 (t = 0.136) and 
Q 2 = 50000 GeV 2 (t = 0.245). We choose to 
work, as in the paper [9], in the fixed flavour 
scheme, nj — 3, with Aq CD — 231 MeV, and 
without taking into account, in the Q 2 evolu- 
tion of a s , quark thresholds. The range ]0, 1] 
has been divided in M steps by M + 1 points: 
xo, x\, xm distributed in such a way that the 
function ln(x) + 2x varies by the same amount at 
any step; this function is slightly different from 
the pure logarithmic distribution commonly used 
in literature [9, 10], but allow, in our case, a more 
uniform distribution of the numerical errors. The 
end points are fixed to be xo = 1 x 10~ 8 and 
xm = 1; however, for a better reading, in the 
Figures 1—4 the x— axis ranges from 10~ 4 to 1. 

First, we fix M = 100. In Figs. 1-2 arc 
reported the evolved singlet and gluon distribu- 
tions, respectively, obtained with n = 3, 6 and 12 
for Q 2 = 200 GeV 2 and Q 2 = 50000 GeV 2 . It 
is worth to note the very fast convergence of the 
series to the solution, as already observed above. 
As a matter of fact, the maximum difference be- 
tween the solutions relative to n = 6 and n = 12 
is 1.4 x 10~ 5 (7.6 x 10~ 4 ) for the singlet, and 
9.3 x 10~ 5 (5.3 x 10~ 3 ) for the gluon distribu- 
tion, in correspondence of Q 2 = 200 GeV 2 (Q 2 = 
50000 GeV 2 ). 

Next we fix n = 12 and Q 2 = 200 GeV 2 . In 
Figs. 3—4 are plotted the approximated evolved 
distributions with M = 25, 50, 100: the maxi- 
mum difference on the common points between 
M = 50 and M = 100 is 4.6 x 10~ 4 for the sin- 
glet and 7.9 x 10 -4 for the gluons. By comparing 
the results in Fig. 3—4 with the corresponding 
Figs. 1-4 in Ref. [9], we observe, besides a good 
numerical agreement of the results, a faster con- 
vergence as the number M of integration points 
increases, as a consequence of our more accurate 
x-integration procedure with respect to the so 
called "brute force" methods. In fact it should 
be observed that reducing from the cubic to the 
linear approximation of f(x) in Eq. (3.1), the 
accuracy £(x,t) (defined in the sequel) becomes 
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about 1 and 3 order of magnitude bigger, respec- 
tively for singlet and gluon. 

To discuss the degree of accuracy of the method 
in [11] were introduced a global accuracy £(x,t) 
defined as the difference between left and right- 
hand side of the Eq. (2.2). The comparison be- 
tween the range of values of £ with the one of 
both sides of Eq. (2.2) represents a very good es- 
timate of the degree of accuracy of the solution. 
In Figs. 5 and 6 are plotted, for n = 12, M = 100 
and Q 2 = 200 GeV 2 both sides of the Eq. (2.2) 
and the corresponding (rescaled) accuracy £(x, t) 
It appears evident that an excellent approxi- 
mation of the solution is obtained. 

Another advantage of our method, once fixed 
the accuracy of the solution, appears to be the 
running time to get each evolution. In fact, the 
simple analytical structure of the evolution ma- 
trix T makes the solution procedure considerably 
fast. As a matter of fact, once given the Split- 
ting Functions and constructed the correspond- 
ing matrices M( fe ) (we have used Mathematica 
[14] to do this), a single evolution, i.e. the multi- 
plication of the T evolution matrix by the initial 
vector, require, for n = 12 and M = 100, about 
6 x 10~ 2 sec on an AlphaServer 1000 using a For- 
tran Code. 

Particularly interesting is the comparison be- 
tween our method and the one presented in [2, 4], 
where an evolution operator is also introduced. 
Firstly we observe that the latter method is based 
on a polynomial expansion of the splitting and 
distribution functions. The expansion is equiva- 
lent to an expansion in power of x. As a conse- 
quence it is affected by problems of convergence 
for x — > 0, due to the branch point in zero of 
the involved functions. This is the source of the 
difficult encountered in the small- a; region, which 
are not present in our approach in which an op- 
timized Newton-Cotes-like quadrature formula is 
employed. 

Second, also in the x-region of convergence, the 
Laguerre polynomial expansion need, for each 
evolution process, the computation of the mo- 

§The integration in the right-hand side has been per- 
formed numerically after an ^-interpolation of the dis- 
crete values obtained with the evolution operator, while 
the left-hand side is worked out by direct derivation of 
Eq. (2.14). 



mcnts of the initial parton distributions with re- 
spect to the polynomials: this procedure requires 
a remarkable amount of CPU-time with respect 
to our approach in which only the evaluation 
of the initial parton distribution in the M grid 
points is needed. 

5. Conclusions 

We have discussed a new algorithm to solve the 
DGLAP evolution equations, in the x-space, which 
appears suitable for a rather large class of cou- 
pled integrodiffcrential equations. 

The method produces a solution which is an- 
alytical in the (Revolution parameter and ap- 
proximate, but rapidly convergent, in the x— space. 
It allows to construct, once for all, an evolution 
operator in matrix form. It depends only on 
the splitting functions appearing in the equations 
and can be rapidly applied to whatever initial 
distribution to furnish the evolved one, requir- 
ing for each evolution only a few hundredth of 
second. 

It is worth to note the reliability of our x- 
integration algorithm, which gets excellent ap- 
proximations on the whole x— range (we use, for 
all the calculations, 10~ 8 < x < 1), also with 
few integration points, resulting in an evolution 
matrix of particularly small dimensions. 

In conclusion, our method, whose numeri- 
cal implementation is straightforward, appears to 
be very fast, very accurate and extremely stable 
with respect to the increasing of convergence pa- 
rameters (i.e. n, the order of the truncated series 
which gives the Evolution Operator, and M, the 
number of integration points). For these reasons 
it represents a powerful tool to analyze the ex- 
perimental data on nucleon structure functions. 
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Figure 1: The initial Singlet distribution (Q — 
4 GeV 2 , solid line) and the evolved ones for n = 
3 (dashed lines), n = 6 (dotted lines) and n = 12 
(solid lines) corresponding at Q 2 = 200 GeV 2 and 
Q 2 = 50000 GeV 2 . We use M = 100. 
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Figure 2: The same in Fig. 1 for Gluons. 
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Figure 3: The initial Singlet distribution (Q 2 = 
4 GeV 2 , solid line) and the evolved one at Q 2 = 
200 GeV 2 with M = 100 (solid line), M — 50 (dotted 
line) and M = 25 (dashed line) with n — 12. 
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Figure 5: For n = 12, M = 100 and Q 2 = 200 GeV 2 
both sides of the Eq. (2.2) are plotted (solid line and 
dotted line), in correspondence of the Singlet distri- 
bution. Dashed line represents £{x) x 10 3 (see text). 
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